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I . INTRODUCTION 


It  is  well  established  that  cross-field  flows  of 

neutral  wind  govern  the  physics  of  plasma  structure  in  the 

high-altitude  nuclear  environment.^  The  description  of 

this  wind  is,  obviously,  a central  issue  in  constructing 

numerical  simulations  of  the  plasma  phenomenology.  In 

many  cases,  limitations  have  arisen  in  previous  research 

and  systems  analysis  efforts  because  of  the  lack  of  a 
, 2 

suitable  routine  to  predict  these  flows.  The  development 
of  a procedure  to  correct  this  deficiency  is  the  subject 
of  this  document.  The  work  is  specifically  directed  at 
constructing  a driver  for  electrostatic  ion  motion  codes 
in  the  late-time  communication  regime.  Oriented  towards 
this  particular  purpose,  the  program  has  a two-fold  objec- 
tive. First,  it  focuses  on  our  desire  to  extend  the 
ability  to  calculate  to  very  late  times  (hours)  and  to 
very  great  distances  (thousands  of  kilometers).  Second, 
it  emphasizes  the  development  of  quick  running,  flexible 
procedures  for  application  in  arbitrary  scenarios. 

The  flow  modelling  is  based  on  the  superposition  of 
acoustic  gravity  wave  (AGW)  modes.  The  heart  of  the 
computational  work  is  a comprehensive,  three-dimensional 
(32x32x32)  AGW  simulator  developed  in  a previous  contractual 
effort.  The  relevant  AGW  mechanics  and  technical  details 
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of  the  code  have  been  thoroughly  described  in  DNA  3420T^ 
and  will  not  be  repeated  here.  Our  purpose  in  this  docu- 
ment is  to  review  the  further  application  of  this  technique 
to  several  new  facets  of  the  modelling  effort  which  have 
been  undertaken  in  the  current  year.  It  does  seem  appro- 
priate at  this  point,  however,  to  quickly  review  the 
rationale  for  selecting  an  AGW  approach  to  the  general 
problem  of  modelling  an  extended  space/time  neutral  wind. 

Following  a nuclear  event,  the  atmosphere  is  set 
into  motion  by  pressure  enhancements  which  arise  either 
directly  at  a burst  or  at  distant  locations  due  to  radiation 
deposition.  For  large  bursts  at  high  altitude,  the  resulting 
motion  is  strongly  influenced  by  the  Earth's  gravity,  both 
directly  and  through  the  stratified  ambient  density.  The 
atmosphere  is  an  acoustic-gravitational  medium  which 
relieves  these  nonequilibrium  pressures  by  sending  out 
pulses  which  set  the  fluid  in  motion.  A natural  represen- 
tation of  the  response  is  in  terms  of  a set  of  acoustic- 
gravity  waves.  In  general,  the  solution  of  the  fully 
nonlinear  fluid  equations  using  the  AGW  approach  is  a 
complex  undertaking.  In  fact,  prior  to  the  advent  of 
Fast  Fourier  Transform  (FFT)  techniques,  one  would  not  choose 
to  go  this  route  except  in  very  special  circumstances.  Thus, 
it  has  become  customary  to  solve  nuclear  event  hydrodynamic 
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[ problems  using  finite  differencing  in  physical  space  of 

the  familiar  fluid  equations.  This  differencing  is 

i - 

' arbitrary  and,  perhaps,  to  a mathematician  "unnatural", 

but  has  been  the  only  real  avenue  open  for  obtaining 
practical  results.  The  introduction  of  the  FFT  has  changed 
the  traditional  odds  markedly.  Today  it  can  be  argued 
that  the  solution  of  fluid  problems  in  appropriate  wave- 
number  space  is  no  more  complex  than  the  older  techniques 
and  stands  to  give  better  answers.  Nonetheless,  we  would 
not  be  prepared  to  make  the  case  that  a complete  AGW  code 
could  improve  on  MRHYDE  or  MICE  without  extensive  and  costly 
development.  It  is  important  to  keep  this  last  statement 
in  mind.  There  is  no  suggestion  in  this  paper  that 
present  general  hydro  codes  are  obsolete  or  ought  to  be 
replaced.  The  discussion  is  concerned  with  augmenting 
the  utility  of  these  codes  by  adding  to  our  capability  in 
an  extended  space/time  regime. 

When  we  shift  our  attention  from  the  early- time/ 
close-in  regime  to  the  late-time/far-distance  regime,  the  ii’ 

equations  that  govern  the  neutral  background  flows  take  on  i 

an  increasing  linear  character.  It  is,  as  we  approach  this  i 

i 

regime,  that  present  hydro  codes  begin  to  experience  serious  I 

difficulties.  It  can  be  argued  that  the  very  characteristics  I 

I 

which  enabled  them  to  do  an  excellent  job  at  the  start  of  | 

i 

I 
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the  calculation  are  now  causing  them  trouble.  The  inclusion 
of  all  the  nonlinear  terms  is  amplifying  noise  to  the 
point  of  unreliability.  The  finite  difference  mesh  forces 
intolerably  small  time  steps.  What  is  clearly  needed  is 
a simpler  reformulation  of  the  problem  which  can  permit 
extension  into  the  new  regime. 

In  the  truly  linear  regime,  such  as  the  ambient 
atmosphere,  there  can  be  little  argument  but  that  flow 
problems  ought  to  be  formulated  in  terms  of  AGW.  The 
technique  has  been  highly  developed  and  is  universally 
used  in  research  concerning  travelling  ionospheric  dis- 
turbances and  other  high-altitude  neutral  phenomena.  The 
popularity  of  the  procedure  is  based  on  the  issue  of 
"natural"  representation.  If  the  medium  can  be  approxi- 
mated as  linear,  then  we  can  express  a disturbance  as  a 
set  of  acoustic  gravity  waves  whose  frequency  and  amplitude 

are  constants . It  becomes  a simple  matter  to  reconstruct 

b 

the  evolving  disturbance  at  any  arbitrary  su^equent  time 
or  place  by  a simple  superposition  of  these  waves.  Utilizing 
the  FFT  to  perform  the  initial  wave  analysis  and  subsequent 
reconstructions,  we  have  a procedure  that  is  quick,  con- 
venient, and  reliable.  In  the  linear  regime,  the  identi- 
fication and  preservation  of  the  constants  of  the  wave 
structures  assures  fidelity  and  absence  of  numerical  noise. 
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The  ability  to  superpose  uncouples  us  from  artificial 
constraints  related  to  grid-size/time-step  requirements. 

It  seems  obvious  to  us  that  a complete  cycle  for 
the  neutral  flow  calculation  ought  to  begin  with  the 
familiar  hydro  code  at  early  time  and  end  with  an  AGW 
code  at  late  time.  Now  if  we  are  to  marry  the  early 
time  and  the  late  time,  there  must  be  a common  meeting 
ground  established  with  a comfortable  overlap.  Having 
commented  on  the  problems  of  nonlinear,  finite  difference 
hydro  codes  in  going  to  late  time,  it  is  appropriate  to 
comment  on  what  happens  to  an  AGW  code  when  it  goes  to 
early  time.  If  we  formulate  our  fluid  dynamics  as  a set 
of  acoustic  gravity  waves  and  monitor  their  behavior  in 
the  highly  nonlinear  regime,  we  will  observe  that  their 
amplitude  and  frequency  are  not  constants,  but  vary  with 
time.  Expressed  more  precisely,  we  will  find  that  the 
evolution  of  a disturbance  does  not  retain  a fixed  modal 
representation  that  can  arbitrarily  be  reconstructed  by 
superposition  of  the  initial  waves.  Instead  we  find  a 
continuous  transfer  of  energy  from  one  set  of  modes  to 
another,  which  alters  the  respective  amplitudes.  Simul- 
taneously, because  the  modes  are  interactive,  the  frequency 
is  shifted  from  its  "natural"  value.  The  calculation  of 
these  amplitude  and  frequency  changes  would  be  the  heart 


of  a fully  coupled  AGW  code.  As  implied  before,  we  have  not 
constructed  such  a code,  but  would  expect  it  to  be  equal 
in  complexity  to  the  traditional  codes. 


.1 
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The  basic  plan  to  marry  a "constant  coefficient" 

AGW  code  to  a general  hydro  code  run  is,  perhaps,  obvious 
at  this  point.  We  need  to  monitor  the  AGW  parameters  in 
the  hydro  code  run  to  determine  when  they  have  settled 
down  and  appear  to  be  approaching  constant  values.  Clearly, 
they  will  never  be  exactly  constant.  It  will  be  necessary 
to  employ  some  engineering  criterion  with  respect  to 
statements  such  as  "sufficiently  constant".  The  next 
section  of  this  review  is  devoted  to  the  subject  of  inter- 
facing an  AGW  solution  to  the  end  point  of  a MICE  hydro 
code  run.^ 

Strictly  from  a physics  veiwpoint,  the  main  task  in 
utilizing  constant  coefficient  AGW  mechanics  for  neutral 
wind  modelling  is  solely  the  question  of  coupling  to  a 
suitable  set  of  initial  conditions.  In  a multiburst  scenario, 
there  is  the  need  to  stop  and  restart  with  new  coefficients 
for  each  new  burst.  However,  this  is  really  just  a compu- 
tational detail  that  was  reasonably  well  explored  in  previous 

3 

work  (DNA  3420T) . Essentially,  if  one  can  couple  to  a 
single  burst,  it  should  be  possible  to  also  couple  to 
repeated  bursts  at  later  tiir«'»s  or  different  locations  with 
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similar  procedures.  If  we  were  using  a Fourier  Integral 
of  infinite  (or  Earth  circumference)  dimension,  there  would 
be  no  additional  computational  task  in  the  development 
program.  As  it  is,  we  employ  Fourier  Transforms  of  distinctly 
finite  spatial  extent  and,  in  particular,  use  fast  Fourier 
procedures.  While  we  have,  from  time  to  time,  looked  into 
alternative  numerical  techniques,  we  have  always  returned 
to  the  position  that  the  use  of  the  FFT  is  an  essential 
ingredient.  It  is  the  FFT  which  permits  extremely  short 
running  time  and  reliable  manipulation  of  massive  blocks 
of  physical  and  Fourier  data. 

The  use  of  FFT  procedures  introduces  a potential 
limitation  on  the  maximum  size  of  our  spatial  grid.  As 
our  objective  is  to  extend  the  flow  solution  to  the  large 
space/time  regime,  it  is  our  desire  to  minimize  any  restric- 
tion of  this  type.  Consequently,  a significant  effort  has 
been  made  to  develop  auxiliary  procedures  to  overcome  this 
problem  with  the  FFT.  We  believe  that  the  current  research 
program  has  been  successful  in  generating  a solution  to 
this  difficulty.  The  third  section  of  this  document  is 
devoted  to  a discussion  of  this  issue. 

We  can  summarize  the  rationale  and  features  of  the 
overall  program  as  follows.  First,  an  extended  space/time 
neutral  wind  model  is  a requirement  of  the  plasma  phenomenology 
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that  is  relevant  to  communication  systems  analysis. 

Second,  a formulation  based  on  AGW  mechanics  is  the  logical 
approach  to  model  construction.  Third,  a basic  AGW  simu- 
lation has  been  built  and  tested  in  earlier  work.  Fourth, 
in  application  there  are  two  development  problems: 

(A)  coupling  to  a hydro  code,  and  (B)  overcoming  any  grid 
limitations  imposed  by  the  use  of  FFT. 
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II.  HYDRO  CODE  COUPLING 


In  order  to  initiate  a legitimate  high-altitude 
nuclear  calculation  using  the  constant  coefficient  AGW 
procedure,  it  has  appeared  desirable  to  start  from  hydro 
code  data  as  indicated  in  the  last  section.  Ultimately, 
for  use  in  broad  ranging  scenario  studies,  it  would  seem 
desirable  to  construct  a family  of  "initializers".  These 
would  constitute  a library  of  simpler  data  sets  that  would 
give  equivalent  flow  properties  for  starting  not  only  an 
extended  wind  model,  but  also  a variety  of  concurrent  late- 
time plasma  modelling  schemes.  Sowle  of  Mission  Research 

Corporation  has  proposed  such  an  activity,  and  we  strongly 

5 

support  the  idea.  In  the  meantime,  the  employment  of 

actual  hydro  code  data  tapes  appears  to  be  the  easiest 

way  to  both  start  AGW  computations  and  to  insert  mode 

corrections  in  simulating  multiburst  flows. 

The  input  and  manipulation  of  hydro  code  data  for 

use  in  extended  space/time  wind  studies  is  conceptually 

straightforward,  but  non-trivial  in  practice.  Fajen  and 

co-workers  at  MRC  have  been  most  cooperative  in  making  the 

library  of  MICE  tapes  available  to  us  and  providing 

4 

instructions  on  reading  and  storing  the  data.  With  their 
assistance,  we  have  become  proficient  in  using  the  data. 

As  a prototype  run  for  use  in  developing  procedures,  we 
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have  focused  on  the  S200  data  with  vertical  magnetic  field. 
This  case  possesses  two  important  features:  1)  it  is  the 
"canonical"  high-altitude  burst,  and  2)  it  is  expected  to 
display  very  large  nonlinear  effects  at  early  time.  In 
understanding  the  behavior  of  this  run,  we  have  assumed  that 
we  were  observing  characteristics  that  would  be  typical 
of  all  such  data  sets. 

The  major  nonlinear  effect  in  high-altitude  bursts 
is  the  rapid  vertical  rise  of  fluid  in  columns  near  the 
fireball  core.  Because  the  neutral  fluid  is  partially 
coupled  to  the  plasma  in  the  core  (which  is  in  turn  tied 
to  the  field) , this  effect  is  maximized  in  a vertical 
field  geometry.  The  altitude  regime  of  primary  interest 
to  us  is  the  space  form  125  km  to  500  km.  The  S200  is 
of  a size  and  location  to  produce  a large  effect  in  this 
range . 

The  first  step  in  building  an  interface  between 
hydro  code  data  and  the  AGW  procedure  is  to  investigate 
the  behavior  of  the  early-time  data  in  terms  of  its 
acoustic  gravity  wave  properties.  We  can  accomplish  this 
task  by  computing  the  actual  mode  amplitude  of  each  AGW 
as  a function  of  time.  Each  mode  is  identified  by  a 
complex  amplitude  which  we  can  write  as 


where  R is  a real  modulus  and  0 is  a real  phase  angle. 
We  can  describe  the  phase  variation  in  time  by  defining 
a real  frequency,  fi. 


6 = Sq  + nt  . (2) 

In  the  linear  regime,  we  would  select  R and  6^  as  initial 
constants.  would  be  a constant  frequency  for  the 
particular  mode  determined  by  a dispersion  relation.  We 
would  determine  the  value  of  A at  any  subsequent  time  by 
merely  consulting  these  constants.  In  the  nonlinear  regime, 
however,  we  can  expect  R to  vary  in  time  and  the  value  of 
6 will  not  be  given  exactly  by  (2).  Thus,  our  specific 
task  is  to  find  the  actual  values  of  R that  MICE  is 
generating  and,  likewise,  the  actual  values  of  0.  With 
respect  to  the  latter,  it  is  convenient  to  remove  the 
natural  phase  rotation  and  consider  only  the  phase  dif- 
ference as  a function  of  time. 

(})  = A0  = 0 - fit  (3) 

The  algebraic  steps  in  constructing  a set  of  complex 
amplitudes.  A,  from  an  arbitrary  data  set  is  straightforward 
but  tedious.  It  has  been  explained  in  detail  in  DNA  3420T^ 
and  will  not  be  repeated  here. 
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As  a general  proposition,  we  can  resolve  the  MICE 
flow  field  into  an  arbitrarily  large  number  of  AGW  modes. 
As  a practical  matter,  it  would  not  make  sense  to  consider 


modes  with  characteristic  dimensions  less  than  the  MICE 
spatial  resolution.  Such  modes  would  prove  to  have 
negligible  amplitude  and  their  consideration  would  be 
wasteful  of  our  computer  storage  capability.  We  have 
found  that  there  is  no  difficulty  in  fitting  the  MICE 
output  if  we  confine  ourselves  to  modes  whose  quarter-cycle 
dimension  is  greater  than  45  km. 

The  other  side  of  the  coin  concerns  the  size  of  the 
maximum  dimension  to  consider  in  the  problem.  To  follow 
the  evolution  of  the  hydrodynamics  for  extended  periods 
in  time,  the  spatial  dimensions  can  become  enormous.  Our 
ultimate  objective  is  to  be  able  to  work  in  a space  that  is 
12,000  km  in  extent.  The  techniques  for  computing  in  such 
a vast  spatial  domain  require  special  development  and  form 
the  subject  of  the  next  section  of  this  document.  For  the 
immediate  purpose  of  investigating  the  MICE  data,  however, 
we  have  not  employed  these  procedures,  nor  are  they  necessary. 
In  the  relatively  early  time  frame  for  which  MICE  computes 
output,  we  find  that  spatial  dimensions  of  3000  km  hori- 
zontal and  1700  km  vertical  are  adequate.  We  therefore 
limit  the  resolution  into  AGW  modes  whose  half-cycle 
dimension  is  less  than  these  values. 
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Bounded  by  the  limits  discussed  above,  our  basic  j 

I 

I 

AGW  simulation  will  construct  amplitude  parameters  for  ] 

i 
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131,072  modes.  The  S200  run  displays  spatial  data  from 
zero  to  about  550  seconds  at  approximately  30  distinct 
time  steps.  Our  task  has  been  to  construct  the  amplitude 
parameters  for  each  mode  at  each  time  step  and  plot  the 
results.  For  each  mode,  we  have  computed  the  time- 
dependent  value  of  R and  (j) . This  result  is,  obviously, 

262,144  separate  curves.  These  have  been  computer  plotted 
and  stored  on  microfiche.  Both  parameters  are  displayed 
in  a common  space  with  banks  of  six  modes  fitted  on  a fiche 
grid,  as  shown  in  Figure  1.  With  this  compression,  we  can 
fit  about  50,000  plots  on  a single  fiche. 

While  the  AGW  data  set  for  the  MICE  S200  run  may 
sound  overwhelmingly  large,  it  is  actually  simpler  to  view 
and  comprehend  than  many  displays  we  have  seen  of  the 
physical  flow-field  data.  Because  all  of  the  curves  are 
essentially  similar,  it  was  possible  to  survey  a vast  array 
simultaneously  and  deduce  the  information  that  we  were 
looking  for. 

Figure  1 is  a convenient  example  to  employ  in  describing 
the  analysis  of  AGW  behavior.  It  happens  to  contain  the 
largest  amplitude  modes  and  thus  was  the  most  important 
single  data  set.  It  displays  the  characteristic  time 
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of  "numerical  noise"  in  MICE,  which  is  not  real  and  should  i 

be  ignored.  In  any  event,  there  is  the  possibility  that  I 

I 

much  of  this  drifting  will  cancel  out  when  the  modes  are  ] 

summed  up  to  obtain  true  physical  representations.  I 

i 

Ideally,  we  would  always  like  to  couple  the  constant  s 

I 

coefficient  AGW  solution  to  a hydro  code  run  at  a point  | 

in  time  when  we  were  absolutely  certain  that  R and  (t  were  i 

I 

stationary.  We  could  then  extend  the  solution  as  far  as  | 

we  wished,  knowing  that  our  numerical  simulation  was  "exact". 

Our  study  of  MICE  S200  has  demonstrated,  however,  that  we 
are  not  likely  to  be  given  such  data.  Undoubtedly,  we 
will  also  have  to  start  when  these  parameters  are  still 
drifting,  the  "quasi-linear  regime".  The  question  will  then 
arise  as  to  the  "legitimacy"  of  coupling  at  this  point. 

Obviously,  we  cannot  provide  a definitive  answer  by  simply 

analyzing  mode  amplitude  data.  The  latter  is  highly  ! 

useful  for  defining  regimes  and  indicating  possibilities. 

It  is  not  clear,  though,  that  it  provides  a simple  quanti-  j 

{ 

tative  assessment  of  the  influence  of  parameter  deviations  | 

on  the  ultimate  physical  solution. 

The  basic  objective  of  the  extended  neutral  wind  1 


program  is  to  provide  cross-field  wind  velocities  and 
neutral  densities  for  use  by  plasma  phenomenology  routines. 
In  the  end,  the  program  will  be  judged  on  the  reliability 
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of  these  results  and  not  on  the  stationarity  of  AGW  mode 
parameters.  It  has  seemed,  therefore,  that  a valuable 
exercise  would  be  to  compare  a constant  coefficient 
simulation  in  the  quasi-linear  regime  to  the  MICE  simu- 
lation. The  time  interval  from  300  to  450  seconds 
appears  to  be  an  appropriate  span  to  consider.  As 
indicated  in  the  discussion  above,  we  believe  that  the 
AGW  mode  analysis  is  telling  us  that  this  is  not  an  un- 
reasonable choice.  Prior  to  300  seconds,  the  S200  appears 
to  show  mode  coupling.  After  450  seconds,  we  observe 
boundary  problems  in  the  MICE  grid  that  are  irrelevant 
to  the  real  physics  that  we  wish  to  simulate. 

In  Figure  2,  we  show  contour  plots  of  the  logarithm 
of  the  density  perturbation  from  the  MRC  MICE  S200  at  300 
seconds.  The  vertical  dimension  is  altitude  in  units 
of  200  km.  The  top  boundary  is  1200  km.  The  "X"  dimension 
of  our  3-D  code  is  equivalent  to  their  radius  in  the  same 
units  as  above.  The  radial  boundary  is  1050  km.  We  are 
going  to  analyze  the  M.RC  data  at  this  time  frame  into  AGW 
modes.  For  two  reasons,  we  do  not  want  their  large  per- 
turbation quantities  up  against  the  boundary  of  our  grid. 

We  desire  to  provide  a reservoir  region  surrounding  the 
data  to  avoid  artificial  steep  gradients  in  the  AGW  reso- 
lution. We  also  want  to  provide  space  for  the  solution 
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Figure  2,  MRC  density  - 300  seconds. 
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to  evolve  without  hitting  the  boundary.  To  accomplish 
this,  we  establish  a vertical  space  from  -50  to  +1750  km 
and  a horizontal  space  from  -1500  to  +1500  km. 

In  Figure  3,  we  show  the  extension  of  MRC  data  into 
this  space  to  provide  a smooth  decay  of  the  perturbed 
quantity.  On  the  same  plot  we  display  the  AGW  fit  to  MRC 
at  300  seconds.  Near  the  edges,  the  AGW  solution  has  been 
additionally  smoothed.  Thus,  it  does  not  fit  exactly  in 
: these  volumes  of  space.  Having  once  fit  the  data  at 

300  seconds,  we  are  going  to  let  the  AGW  code  run  with 

constant  coefficients  and  no  further  interaction  with  the  j 

MRC  data.  We  can  run  the  comparison  out  to  450  seconds. 

Figure  4 shows  the  comparison  at  330  seconds,  while  \ 

in  Figure  5 we  jump  out  to  450  seconds.  Figures  6,  7,  and 
8 provide  the  same  display  for  the  radial  velocity.  The 
density  and  cross-field  velocity  are  the  two  parameters 
fed  to  the  plasma  electrostatic  integration.  We  could, 

j 

of  course,  generate  any  other  flow  parameter  if  there  was 
a need. 

It  is  important  to  emphasize  that  the  figures  show 
both  the  battle  space  and  the  reservoir.  The  reservoir  is 
garbage  for  purposes  of  this  comparison.  One  should  review 
Figure  2 to  observe  the  actual  MICE  grid  space  which  would 
be  equivalent  to  the  battle  space.  We  believe  that  the 
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comparison  in  the  battle  space  region  is  quite  encouraging. 

We  do  note  that  there  are  problems  in  getting  a 
good  comparison  along  the  central  axis  core.  The  AGW 
solution  has  a tendency  to  overemphasize  the  central 
ballistic  rise  of  fluid.  (It  is  well  known  that  "linear" 
procedures  overestimate  effects  in  the  nonlinear  regime 
because  of  the  lack  of  the  convective  derivative,  which 
typically  provides  a saturation  or  damping  influence.) 

We  find  the  density  contours  a little  higher  in  this 
space,  while  the  radial  velocity  contours  are  distorted 
by  the  upsweep  of  fluid.  We  have  made  no  attempt  to 
"tune"  our  simulation  to  fit  the  MICE  data.  We  could 
probably  juggle  the  gas  properties  or  ambient  atmosphere 
to  correct  some  of  the  central  core  discrepancy.  We  do 
not  believe  that  this  would  represent  the  proper  spirit 
of  the  simulation,  however. 

Away  from  the  center  line,  the  flow  properties  are 
reasonably  well  simulated.  This  is,  after  all,  the  major 
region  for  striation  development  and  represents,  by  far, 
the  greatest  volume  (which  may  not  be  obvious  in  a 2-D 
plot)  . 

Having  coupled  to  the  MICE  run  and  computed  the  AGW 
parameters,  we  can,  of  course,  run  the  solution  out  in  time. 
Figure  9-14  show  the  extension  to  2400  seconds.  In  this 
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Figure  14.  AGW  radial  velocity  - 2400  seconds 
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case,  we  have  sliced  off  the  reservoir  in  the  display  to 
produce  an  easier-to-read  result.  We  do  not  attempt  to 


run  beyond  this  point  because  of  boundary  problems  of  our 
own.  Particularly,  in  the  velocity  plot  at  2400  seconds, 
we  observe  what  appears  to  be  artificial  noise  that  we 
believe  is  coming  from  boundary  reflections. 

In  reviewing  the  work  described  in  this  section,  we 

can  identify  two  areas  for  further  development.  First,  it 

would  be  desirable  to  start  from  hydro  code  data  that  does 

not  have  large  disturbances  up  against  its  own  grid 

boundaries.  This  would  eliminate  the  need  for  the  somewhat 

wasteful  "reservoir  region".  For  the  specific  case  of 

S200,  AFWL  is  preparing  such  data  for  us  using  an  extended 

6 

grid  and  the  HULL  code.  As  a general  proposition,  we 
probably  need  to  pre-process  hydro  code  data  so  that  it  is 
in  better  format  for  our  application.  As  a long-term 
objective,  we  look  forward  to  the  development  of  convenient 
"initializers",  as  mentioned  elsewhere. 

A second  area  of  development  is  clearly  the  problem 
of  "extended  boundaries".  We  do  not  want  to  give  up  our 
45  km  resolution,  but  on  the  other  hand,  we  do  not  want 
to  increase  the  mode  number  above  its  present  value  of 


131,072.  The  curious  reader  may  inquire  as  to  why  we  need 
so  many  modes  in  the  first  place.  Physically,  we  can 


observe  that  there  may  only  be  several  hundred  "key  modes" 
needed  for  an  excellent  fit  to  a particular  data  set. 

Why,  then,  do  we  insist  on  an  absolutely  "complete  set" 
of  all  possible  modes?  The  answer  is  two-fold.  In  the  one 
case,  we  want  to  retain  the  complete  set  of  acoustic 
gravity  wave  physics  in  order  that  we  can  treat  arbitrary 
initial  conditions.  We  will  not  a priori  know  the  regimes 
dominated  by  pressure  waves  (compression  or  rarefaction) 
or  by  in-running  or  out-running  pure  gravity  waves.  In 
the  second  case,  we  wish  to  retain  a total  spatial  resolution 
in  order  to  effectively  utilize  the  FFT.  We  have  explored 
the  concept  of  incomplete  Fourier  sets,  but  these  require 
giving  up  the  FFT.  It  is  important  to  remember  that  it 
is  the  FFT  which  permits  us  to  time  step  this  vast  mode 
array  in  a small  fraction  of  a second  on  the  7600.  We 
believe  that  there  is  an  alternative  solution  to  the  ex- 
tended boundary  question.  This  forms  the  subject  of  the 


next  section. 


III.  EXTENDED  BOUNDARIES 


In  previous  sections,  we  have  indicated  that  the 
use  of  finite  Fourier  Transforms  imposes  distinct  boundaries 
on  the  AGW  modelling  effort.  We  desire  to  extend  these 
boundaries  to  very  large  values.  Specifically,  there  are 
two  related  reasons  for  wishing  to  extend  the  overall 
dimensions  of  the  grid.  The  first  is  the  problem  of  ri<mning 
to  very  late  time  in  the  principal  battle  space  itaelf. 

If  we  do  nothing,  there  will  arise  an  apparent  "reflection" 
of  waves  when  the  disturbances  reach  the  edge.  Cremated 
by  the  periodicity  in  the  finite  Fourier  Transforms,  this 
effect  is  unphysical  and  must  be  removed.  While  various 
"fixes"  can  and  have  been  imposed  in  the  past,  they  in- 
evitably introduce  an  element  of  noise  and  uncertainty 
which  grow  with  time.  The  ideal  solution  is,  of  course, 
sufficiently  extended  boundaries  to  avoid  the  ^oblem  on 
Lime  scales  of  interest.  The  second  reason  for  wishing 
for  more  space  is  the  simple  desire  to  generate  flow  infor- 
mation at  more  remote  locations.  Our  present  data  base 
shows  that  disturbed  wind  fields  (of  interest  to  striation 
analysis)  extend  to  the  edge  of  present  grids — and  thus, 
logically,  beyond. 

Obviously,  the  simple  answer  to  the  grid  problem  is 
to  merely  increase  its  size.  Were  there  no  limits  to 
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computer  storage  this  would,  of  course,  be  done.  The  fact 
is,  though,  that  we  already  exploit  this  parameter  about 
as  far  as  practicable.  We  do  not  wish  to  intrude  on  a 
major  feature  of  the  scheme  which  empahsizes  fast  running 
with  no  off-line  interactive  storage.  Merely  "stretching" 
the  grid,  on  the  other  hand,  would  run  the  risk  of  de- 
creasing the  resolution  below  acceptable  limits. 

It  occurred  to  us  some  time  ago  that  the  inherent 
characteristics  of  the  FFT-AGW  approach  ought  to  allow  a 
sequential,  flexible  data  handling  that  could  get  around 
the  problem  of  large  unit  storage.  It  was  not,  however, 
unitl  the  problem  was  discussed  with  L.A.  Wittwer  that  a 
viable  scheme  was  organized.^  He  suggested  the  concept 
of  analytic  interpolation  of  intermediate  modes  with 
storage  only  of  particular  reference  modes.  Instead  of  a 
straightforward  total  FFT  of  a "big  grid" , he  indicated 
that  it  ought  to  be  possible  to  break  the  problem  into 
pieces  of  manageable  size  by  reorganizing  the  finite 
Fourier  sums.  Acting  on  this  intuitively  attractive 
possibility,  we  have  translated  the  concept  into  a 
rational  mathematical  and  numerical  routine.  An  outline 
of  the  details  is  presented  as  an  appendix  to  this 
document.  We  believe  that  the  concept  which  underlies  this 
approach  will  have  applicability  considerably  beyond  the 
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particular  use  that  we  need  for  the  AGVJ  modelling  problem. 

We  have  tested  the  concept  for  a relatively  simple 
one-dimensional  gravity  wave  example.  The  results  of 
this  test  will  be  described  in  the  figures  which  follow. 

We  have  begun  the  task  of  incorporating  the  scheme  into 
the  full  3-D  AGW  simulator  (the  ACCOU  code) , but  the 
completion  of  this  work  lies  in  the  future. 

In  Figure  15,  consider  a simple  gravity  wave  in 
a 32-grid  space.  As  an  example,  consider  each  horizontal 
tick  mark  to  represent  50  km.  Figure  16  shows  the  con- 
ventional evolution  of  the  wave  in  the  1600  km  total  space. 

Each  vertical  display  represents  the  wave  amplitude  at 
subsequent  time  steps.  As  an  example,  consider  each 
vertical  display  to  be  100  seconds  later  than  the  display 
above.  For  the  1 km/sec  propagation  speed  which  characterizes  < 

this  wave  packet,  we  get  an  apparent  reflection  at  the  J 

boundary  in  800  seconds.  By  the  time  we  have  progressed  | 

to  the  bottom  (2400  seconds) , the  space  is  completely  J 

muddied  up  and  the  physical  description  bears  no  relation-  ] 

ship  to  a real  unbounded  atmosphere. 

Figure  17  is  the  same  battle  space,  but  now  the 
effective  boundaries  have  been  extended  twice  as  far, 
equivalent  to  a 3200  km  total  space.  We  display,  however, 
only  the  original  1600  km  space.  Now  the  wave  has  twice  ! 


as  far  to  go  and,  considering  the  return  time  also,  will 
take  almost  four  times  as  long  to  be  strongly  reflected 
back  to  the  reference  grid.  Figure  18  continues  the 
scenario  out  to  4800  seconds  and  by  this  time,  we  see 
that  twice  as  big  is  not  enough. 

In  Figure  19,  we  now  play  the  same  game  in  grand 
style.  The  effective  boundaries  are  now  made  8 times 
farther  away,  equivalent  to  a 12,800  km  space.  Nothing 
comes  back  in  2400  seconds,  and  Figure  20  shows  that 
nothing  comes  back  in  4800  seconds.  In  each  of  these 
cases,  we  have  only  taken  the  trouble  to  return  the 
original  1600  km  space  with  the  same  50  km  resolution. 

In  Figure  21,  we  show  the  grand  scenario  for  the 
8X  example.  We  go  for  broke  and  return  the  entire  256 
space  using  8 sequential  32-space  FFT's.  The  total  space 
displayed  is  now  the  full  12,800  km  and  we,  of  course, 
see  why  the  reflection  problem  is  non-existent.  We  also 
see  the  flexible  nature  of  the  scheme  which  can  return 
the  big  space  at  late  time  when  it  may  be  needed. 
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Figure  20.  8x32  grid  - time  sequence-2. 
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IV.  FUTURE  DIRECTIONS 


The  most  important  immediate  task  in  AGW  model 
development  is  the  continuation  of  the  effort  to  incorp9j-ate 
the  extended  boundary  procedure  into  a full  3-D  realization. 

On  a longer  term  basis,  the  question  of  efficient 
coupling  to  general  hydro  code  data  will  return.  In 
particular,  we  have  not  attempted  to  interface  hydro  code 
inputs  in  multiburst  scenarios.  Earlier  work  did,  how- 
ever, display  the  capability  of  stopping  and  restarting 
with  an  artificial  set  of  hydro  data.  We  would  like  to 

explore  the  option  of  starting  from  a radiation  code 

2 

such  as  the  ROSCOE  deposition  routine.  This  capability 
would  provide  a wide  range  of  easily  accessed  initial 
conditions  for  parametric  studies.  Ultimately,  we  would 
hope  to  work  closely  with  any  effort  to  develop  a library 
of  flow  "initializers". 

As  the  AGW  model  begins  to  produce  usable  output, 
we  will  need  to  explore  the  interface  requirements  with  a 
body  of  plasma  phenomenology  routines  in  the  community. 
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APPENDIX  - THE  WITTWER-CHU  INTERPOLATION  SCHEME 
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For  simplicity  in  discussion  and  notation,  let  us 
describe  the  interpolation  scheme  in  one  dimension.  The 
extension  to  three  dimensions  is  straightforward  but 
tedious.  Consider  a finite  Fourier  Transform  of  a function, 
f,  which  is  defined  over  some  interval,  L,  and  is  presented 
as  a spatial  quantity  in  the  variable,  x 


f(Xs) 


= 1: 


(k^) 


ik^x 
t s 


(Al) 


The  "s"  is  an  index  for  the  space  grid  of  dimension  Ax  and 
the  "t"  is  an  index  for  the  Fourier  grid  of  dimension  Ak. 


Xg  = s Ax  (A2) 

k^  = t Ak  (A3) 

Ax  = L/N  (A4) 

Ak  = 2Tr/L  (A5) 


1 
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The  resolution  in  physical  space  is  Ax  and  the 
solution  is  periodic  (and  meaningless)  outside  of  the 
interval  L.  The  corresponding  maximum  wavenumber  in 


Fourier  space  is  given  by 


NAk 


and  all  higher  wavenumber  modes  are  assumed  not  to  exist. 


If  f is  a smooth  physical  variable  that  we  have 


properly  resolved  on  a grid  of  dimension  Ax,  the  resolu- 


tion in  Fourier  space  on  dimension  Ak  will  be  satisfactory 


and  the  function  f will  be  smooth  in  that  space. 


Because  we  have  a smooth,  well  resolved  function. 


we  would  have  no  trouble  interpolating  for  values  of  f 


on  dimensions  smaller  than  ax  between  grid  points  in 


physical  space.  Likewise,  for  similar  reasons,  we  will 


have  no  trouble  interpolating  for  values  of  f on  dimensions 


smaller  than  Ak  between  grid  points  in  Fourier  space, 


The  information  obtained  by  interpolating  f in 


Fourier  space  is  exactly  the  data  required  to  extend  the 


boundary.  We  can,  in  fact,  extend  the  boundary  as  far  as 


we  wish,  if  we  are  willing  to  increase  our  computing  time 


and  allow  for  some  additional  temporary  storage, 
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Let  us  assume  that  we  want  to  extend  the  boundary 
(or  periodicity)  by  M times,  but,  of  course,  wish  to  keep 
Ax  fixed.  The  new  grand  scale  dimension  becomes 


L'  = ML 


(A8) 


and  the  corresponding  increment  in  Ak  becomes 


6k  E (Ak) ' = = i (Ak) 


(A9) 


We  can  now  rewrite  our  Fourier  sum  in  (Al)  as 


N-1  M-1  ,,  , „ , 

^ ^^^t  e (AlO) 


t=0  J=0 


or  with  a redefinition  as 


M-1  r N-1  • 1 ■ 

Z "'t'  ® 

£=o  L t=o  -I 


e6k 


(All) 


where 


g,,=o  <^t^  = 


(A12) 


and 


^ f + «.6k) 


(A13) 


■1 
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We  note  that  the  inner  sum  in  (All)  can  be  expressed 


as  a new,  artificial  physical  space  function,  g , 


N-1 


ix  k. 


Z/s  A ^ , 

(k^)  e 


(A14) 


t=0 


Thence  we  can  express  the  grand  scale  function  in 
(All)  as 


M-1 


■n  = Z ^^n^  ® 


ix  {,6k 
n 


(A15) 


^=0 


At  this  point,  we  have  assembled  all  of  the  basic 
tools  that  are  required.  Let  us  walk  through  a computational 
step  to  see  what  actually  goes  on.  Equation  (A12)  represents 
the  stored  array  of  Fourier  data  that  is  always  available 
and  which  we  would  have  used  directly  had  we  not  developed 
the  interpolation  scheme.  Equation  (A13)  represents 
Fourier  data  that  is  always  simply  derivable  from  (A12) 
by  interpolation.  It  is  never  permanently  stored  and  is 
created  only  when  required. 

Equation  (A14)  represents  an  FFT  operation  to  create 


the  g.  array.  Note  that  it  represents  an  N-dimensional 


operation.  The  further  step  to  equation  (A15)  which 
creates  the  real  physical  variable  f is  a subsequent  FFT 
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operation  on  the  array.  Note  that  it  represents  an 
M-dimensional  operation. 


We  would  term  the  above  sequence  of  steps  the  basic 
operational  mode.  We  have  reconstructed  our  basic  physical 
variable,  f,  in  the  original  space  of  dimension,  L.  If 
we  were  reconstructing  f from  our  array  of  f as  a flow 
variable  at  "early  time",  we  would  presumably  see  no  dif- 
ference from  that  of  a conventional  solution.  We  would 
have  done  more  work,  however.  Instead  of  a single  N- 
dimensional  FFT,  we  would  have  required  M number  of  N- 
dimensional  FFT's.  We  would  also  have  had  to  temporarily 
store  twice  as  much  data,  although  the  data  required  in  the 
machine  at  any  moment  to  perform  an  FFT  was  identical. 

If,  on  the  other  hand,  we  were  reconstructing  f 

/V 

from  the  f array  at  "late  time",  we  would  see  a dramatic 
difference  in  the  two  cases.  The  conventional  case  would 
have  suffered  boundary  "reflections"  and  become  meaningless. 
The  new  case  would  have  proceeded  beautifully  just  as  though 
the  boundary  associated  with  L did  not  exist.  It  would  be 
functioning  in  a space  equivalent  to  ML  in  extent.  That  is, 
we  would  only  "see"  the  data  in  the  interval,  L,  but  the 
solution  would  "exist"  throughout  ML. 


Now  suppose  we  had  tried  to  accomplish  the  above 
solution  by  conventional  means.  We  would  first  have  had  to 


I 


I 


generate  and  store  a permanent  f array  that  was  M times 
larger  than  before.  Second,  we  would  have  had  to  bring 
this  data  into  position  all  at  once  to  perform  an  MN- 
dimensional  FFT. 

As  indicated  in  an  earlier  part  of  the  discussion, 
the  effective  removal  of  boundaries  in  the  original  battle 
space  interval,  L,  is  not  the  only  purpose  of  the  new 
procedures.  One  of  our  objectives  is  to  obtain  flow  infor- 
mation at  distant  locations  at  late  time.  To  accomplish 
this  task,  we  employ  the  same  tools,  but  require  one 
additional  step  of  analysis.  A few  paragraphs  back,  when 
we  came  to  Equation  (AlO) , the  careful  observer  may  have 
noted  an  apparent  change  in  notation  without  any  mention 
in  the  text.  The  quantities  f^  and  appeared  in  place 
of  the  expected  f and  x^ . In  the  basic  mode  of  operation, 
these  are,  in  fact,  identical  quantities.  The  new  indices 
take  on  meaning  when  we  consider  a reconstruction  in  the 
extended  domain. 

Let  us  assume  that  we  are  interested  in  reconstructing 
the  flow  in  some  region  of  our  extended  space.  We  desire 
the  information  over  some  interval  of  dimension  L,  but  not 
our  original  interval.  We  define  this  space  as  the  "n^^" 
interval  (one  of  a total  of  M) . Our  basic  physical  space 
grid,  x,  can  be  indicated  in  the  new  interval  by  a simple 
translation 


I 
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(A16) 


Our  previous  expressions,  such  as  (A14) , are 
periodic  in  L,  thus 

We,  thus,  can  return  and  use  (A15)  just  as  before, 
but  now  we  employ  the  new  "x^"  in  place  of  the  as  before. 
Mathematically,  we  create  a phase  change  by  shifting  to  the 
extended  value  for  "x"  and  then  perform  the  same  FFT 
operation.  This  step  enables  us  to  recover  the  real 
physical  variable  f in  an  arbitrary  interval  "n"  which 
was  the  purpose  of  the  notation  f^. 

If  we  are  ambitious,  we  can  reproduce  the  entire  ML 
space.  We  do  this,  however,  by  reconstructing  each  part 
(piece  by  piece)  in  sequence.  This  permits  us  to  develop 
a massive  data  bank,  if  we  desire,  without  increasing  our 
storage  requirement  on  the  core  of  the  computer.  Obviously, 
if  we  want  to  keep  all  of  this  data  for  subsequent  operations, 
it  must  be  stored  outside. 

It  is  by  use  of  the  Wittwer-Chu  procedures  that  we 
hope  to  maintain  50  km  resolution  at  one  point  in  the 
simulation  and  still  consider  an  unbounded  space  in  excess 
of  12,000  km  dimensions. 
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